
/// MLSDEFORMER source code
/// by Alvaro Cuno, LCG-COPPE-UFRJ 
/// bugs? email-me: alvaroecp@gmail.com
/// January, 2007

#ifndef __DEFORMATION_H__
#define __DEFORMATION_H__

#include "../typedefs.h"
#include "constraints.h"

using namespace std;

/// Given two control points sets, this class 
/// applies an optimal deformation using rigid-body transformations
/// on a set of points in a MLS-sense
class Deformation {
   
public:
   Deformation();
   ~Deformation();
   
   /// Sets the data used to perform the MLS deformation
   ///@param ipos pointer to initial positions
   ///@param n number of vertex positions
   ///@param icpos pointer to control points 
   void setting(const vector<Point3> *ipos, unsigned n, const vector<Point3> *icpos);
   
   /// Sets the data used to perform the MLS deformation
   ///@param ipos pointer to initial positions
   ///@param n number of vertex positions
   ///@param icpos pointer to control points
   ///@param W  pointer to table of weights
   void setting(const vector<Point3> *ipos, unsigned n, const vector<Point3> *icpos, Number *W);
   
   /// Precomputes the necessary tables for performing the MLS deformation
   void initialize();
   
   /// Executes the MLS-deformation induced by the deformed control points positions
   ///@param dcpos pointer to current control points (deformed positions)
   ///@param p pointer to current vertex positions (deformed positions)
   void execute_deformation(const vector<Point3> *dcpos, Number *p);
   
   /// Applies the transformations on the precomputed values
   void rotate(const Point3 &o, const Vector3 &u, Number angle);
   void translate(const Vector3 t);
   
   void clear();
   
   /// The weight function
   inline Number omega(Number d);
   
protected:
   /// Auxiliar functions for the deformation algorithm
   Number computeLargestRootOfQuarticEquation();
   void compute_u();
   void computeRotation(const Number v[]);
   void computeCorrelationMatrix_points(unsigned i);
   void compute_pstar_points(unsigned i);
   void compute_qstar_points(unsigned i);
   
   /// Computes the weights table using the Euclidian metric
   void compute_WeightsTable();
   /// Computes auxiliar values: pstar list and vhat list
   void compute_pstarList_vhatList();
   /// Computes the table of phats
   void compute_phatTable();
   /// returns the weight between a vertice and a control point
   inline Number weight(unsigned ivertice, unsigned ictrlpoint) {
      return WeightsTable[ictrlpoint*texSizetexSize + ivertice];
   }
   
protected:
   /// AUXILIAR VARIABLES
   
   /// These variables are valid for the current vertex position computation
   Number                   traceM;     //< trace of M
   Number                   detM;       //< determinant of M
   Number                   eigenvalue; //< eigenvalue
   Number                   pstar[3], qstar[3];
   Number                   Wt, W;      //< weights
   Number                   qjhat[3];
   Number                   u[3];     //< u value 
   Number                   vxu[3], rv[3], factor;     //< a factor value
   Number                   Mf0[3], Mf1[3], Mf2[3]; //< MQ = Sum w phat x qhat   correlation matrix
   Number                   Nf0[3], Nf1[3], Nf2[3]; //< N = MQ + transpose(MQ)
   Number                   V[3];                   //< V = Sum w phat dot qhat
   Number                   WEXPONENT;  //< the exponent of omega function
   
public:
   /// AUXILIAR TABLES
   vector<Vector3>          pstarList;       // 
   vector<Vector3>          qstarList;       // 
protected:
public:
   vector<Vector3>          vhatList;        //< vhat list: v - pstar
   vector<Number>           wtl;             //< w total list: one for each vertex
   
   Number                  *phatTable;          //< phat list: It is used to compute M = w qhat phat. It will be send to GPU
   Number                  *WeightsTable;       //< table of weights
   const vector<Point3>    *i_vpositions;       //< initial positions of vertices that will be deformed
   const vector<Point3>    *i_constraintsPos;   //< deformed positions of contraints
   const vector<Point3>    *d_constraintsPos;   //< deformed positions of contraints
   
   unsigned                 npositions;         //< the number of vertices
   unsigned                 nconstraints;       //< number of control points
   
   unsigned                 texSizetexSize;     //< texture size x texture size
   
protected:
   
   inline Number sgn(Number v) const {
      if (v<0) return -1;
      else return 1;
   }
   
   /// Returns one cubic root of a cubic equation x^3+ax^2+bx+c: Using xxxx
   Number getCubitRoot_X(const Number a, const Number b, const Number c);
   /// Returns one cubic root of a cubic equation x^3+ax^2+bx+c: Using yyyy
   Number getCubitRoot_Y(const Number a, const Number b, const Number c);
   /// Returns one cubic root of a cubic equation x^3+ax^2+bx+c: Using zzzz
   Number getCubitRoot_Z(const Number a, const Number b, const Number c);
   
   /// Constants used for computing cubic roots
   static const Number _1DIV3 = 0.333333; // 1/3
   static const Number _1DIV9 = 0.111111; // 1/9
   static const Number _1DIV54 = 0.0185185; // 1/54
   static const Number SQRT3= 1.73205080; // sqrt(3.0);
   static const Number EPSILON = 1e-5;
};

#endif
